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Abstract- The Iterative Fresnel Integrals Method (IFIM) has been 
applied for the simulation and generation of the complete near- 
field Fresnel diffraction images created by triple apertures for 
the first time. The simulation can be performed in any PC using 
an included MATLAB program. Necessary formalism has been 
derived and simulation algorithm has been developed for this 
application. An interesting combination of interference effects 
with Fresnel diffraction has been observed in the simulation 
images. Transition to the expected Fraunhofer diffraction 
pattern from Fresnel diffraction for triple apertures is also 
observed by the simulations. The program can serve as a useful 
tool to study the complex phenomenon of Fresnel diffraction 
from triple apertures. 
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I. INTRODUCTION 

In Optics, diffraction is a subject of recurring importance, 
both from a theoretical and experimental point of view [1-3]. 
Diffraction phenomena can be grossly classified into two 
types, namely: Fresnel and Fraunhofer. Of these two, the 
Fresnel or near-field diffraction is relatively more complicated 
than the Fraunhofer, or far-field diffraction. No exact 
analytical solution of Fresnel diffraction can be found even in 
the simplest cases. For example, even for a simple rectangular 
or circular aperture, no analytical solution can be obtained for 
Fresnel diffraction. In this case, a solution can be obtained in 
terms of non-analytical integrals known as Fresnel integrals 
[2,3]. Even then, visualizing the complicated Fresnel 
diffraction can be a difficult task. Usually some visualization 
tools, such as Cornu Spiral, are necessary, and even if it's 
used, the complicated diffraction pattern observed is not easy 
to explain and interpret [2]. 

We previously introduced a new method of calculation of 
the complete Fresnel diffraction pattern from rectangular- 
shaped apertures by the Iterative Fresnel Integrals Method 
(IFIM). The technique has been applied to a single rectangular 
aperture [4], double apertures and its derivatives [5] and a 
tilted aperture [6]. Instead of the usual methods of calculation 
of Fresnel diffraction from apertures, which generally uses 
FFT-based algorithms [7-9], the IFIM method uses repeated 
calculation of Fresnel cosine and sine integrals and virtual 
displacements of the aperture [4]. The method provides a very 
intuitive method of calculation of the diffraction pattern, and 
two-dimensional images can be directly generated by the 
method in any arbitrary experimental configurations, such as 
for any given aperture size, illumination wavelength, and 
aperture-screen distances. The algorithm was implemented in 
MATLAB, and the codes can be executed in any PC or laptop 



in a reasonable amount of time, with execution times less than 
a minute in most cases. Implementation of the codes in other 
languages, such as Mathematica or MathCad is also possible. 

In this paper, we extend the IFIM technique to a triple 
aperture. We first derive the basic equation to describe the 
electric field or intensity due to the triple aperture, and then 
formulate the detailed algorithm for the computation process. 
The algorithm is then implemented in MATLAB . We present 
output images from the MATLAB program for typical cases 
of the triple aperture problem. Finally^ discussions are made 
and conclusions are drawn in the final sections. 

n. Theoretical Considerations 
To understand the underlying theory of Fresnel diffraction 
from a triple aperture, let us first consider the case of a single 
aperture. Light of wavelength X emitted from a point source S 
is diffracted by a rectangular aperture of dimension aXb 
located at a distance p from it (Fig.l). The diffracted light is 
observed on the observation plane (or screen) placed a 
distance q away. For convenience, the coordinate systems on 
the aperture and on the image planes are chosen to be centered 
on the optical axis passing through the center of the aperture 
and normal to it, and are denoted by (y, z) and (Y, Z) axes, 
respectively. The Huygens-Fresnel principle is then invoked 
to calculate the total electric field at any given point of the 
image plane (Y,Z) by summing up the contributions, taking 
into account both amplitude and phase, of all the elementary 
wavelets emitted by different area elements inside the clear 
rectangular aperture. 




Fig. 1 Basic configuration of Fresnel diffraction from a rectangular aperture 

The contribution to the complex electric field at the point 
P (located at the origin of the image plane or screen) due to 
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waves emitted by the surface element dS located on the 
aperture is given by [2, 3] 

dE= £ ° K ^\ xp[j{k(p + q) - a> t]dS (1) 
pqA 

where ^is the electric field of the source, k is the 
wavenumber (2n/X) and £jf^M/2(l+cos c|>) is the obliquity 
factor [2, 6], where </> is the angle between the normal to the 
surface element dS and the direction of wave propagation. 
The obliquity factor K((j)) takes into account the decrease of 
intensity for secondary Huygens wavelets emitted from dS in 
an oblique direction. In our case, since q is much larger than 
the dimensions of the aperture, 0is quite small, and K((/)) is 
nearly unity. 

The net electric field at P due to the entire tilted aperture is 
obtained by summing up the contributions dE over the entire 
rectangular aperture, 

E P = I ^Y-exp[j{k(p + q)-co t]dS. (2) 

aperture P ^ 

By integrating the contributions of the Huygens wavelets 
over the area of the entire aperture, it is shown in textbooks [2, 
3] that the total electric field at P is given by 

E p = ^[C(u) + /S(k)]£[C(v) + iS(y)]% (3) 

where E u is the unobstructed electric field at P (i.e. the electric 
field that would have existed if the aperture were absent), and 
C(u) and S(u) are the Fresnel cosine and sine integrals, 
defined by, 



w 

C(w)=j cosOw a I2)dw } and 

o 

w 

S{w) = \sm(nw a l2)dw\ 



(4) 



Here w represents either of the two dimensionless variables u 
or v, being defined as, 



I 2(p + q ) 
X p q 



V A Po% 



The variables u and v are obviously proportional to Cartesian 
coordinates y and z. Specifically, the limits u h u 2 v ly v 2 
correspond to the values of y b y 2 , Zi, z 2 , respectively. 

The intensity at P is given by the square of E p , appearing 
in Eq.(3) i.e., by, 

/„ =±{[C{u 2 ) - C(u t )f + [S(u 2 ) - S( Ul )f } x 

{[C(v 2 )-C(v 1 )f+[5(v 2 )-5(v 1 )f}. (6) 



of the central aperture are then located at y-a/2 and y=-a/2 
respectively, and the edges of the right aperture is located at 
y=a/2+b and y=3a/2+b respectively. Similarly, the edges of 
the left aperture are located at y=-a/2-b and y=-3a/2-b. In this 
notation, the aperture width is a and center-to-center aperture 
separation is (a+b). The upper and lower edges of the 
apertures are located at z=c and z=-c respectively. 

Under illumination from the source, the total electric field 
at the center P of the image plane consists of three 
contributions: one from right aperture, one from the central 
aperture, and the one from the left aperture. The total electric 
field contribution from the right aperture is given, in analogy 
toEq. (3), by 

£ ffl 4[C(tt) + /S(tt)];[C(v) + tf(v)]; (7) 

Here, the limits uj and u 2 are the values of the dimensionless 
variable u corresponding to two edges of the right aperture, i.e. 
for yj=a/2+b and y 2 =3a/2+b, respectively. Similarly, the 
limits vj and v 2 are the values of the dimensionless variable v 
corresponding to the lower and the upper edges of this 
aperture i.e. for Zi=-c and z 2 =c respectively. 

The total electric field contribution by the central aperture 
is similarly given by 

e pc =4[c(«) + iS( U )] nav) + is( V )] (8) 

where u 3 and u 4 are the values of u corresponding to two 
edges of the central aperture, i.e. for y 3 =-?J2 and y 4 =a/2 
respectively. The values vj and v 2 are the same as in Eq. (3) or 
Eq.(7). 

The total electric field contribution by the left aperture is 
similarly given by 

E p =^[C(u) + iS(u)] u ;[C(v) + iS(v)];» (9) 

where u 5 and u 6 are the values of u corresponding to two 
edges of the central aperture, i.e. for y 5 = -(3a/2+b) and y 6 =- 
(a/2+b), respectively. The values v 7 and v 2 are the same as in 
Eq. (7) and Eq. (8). 

The total electric field at P contributed by the triple 
aperture is given by the arithmetic sum of the three complex 
amplitudes from the left, centre and right apertures; i.e. by, 

E = E PR +E PC +E PL 

=^{[C(u) + jS(u)]2+[C(u) + jS(u)]:i+ 

[C(u) + jS(u)]:i } x [C(v) + jS(v)] I (10) 

The intensity at P is proportional to (EE*), i.e., 

I p=\[{C(u 2 ) + C(u 4 ) + C(u 6 )- C(u x )- C(u 3 )- C(u 5 )} 2 
+ {S(u 2 ) + S(u 4 ) + S(u 6 )- Siu^- S(u 3 )- S(u 5 )} 2 ] x 



[{C(v 2 )-C(vJ} 2 +{S(v 2 )-S(vJ} 2 }. 



(ID 



Here, I is the unobstructed intensity corresponding to Here Io denotes the intensity of the unobstructed wave, i.e. Io= 
£ w (Io=E u )• E u . 



Let us now consider the much more complicated case of a 
triple aperture. As before in the single aperture case, we 
assume that the aperture is centered on the yz coordinate 
system, i.e., the origin of the coordinate system O is located at 
the exact center of the triple aperture (Fig. 2). The two edges 



If b is set to zero, then we have a single aperture of width 
3a. In that case, both u 4= ui and u 6= u 3 , and, hence only C(u 2 ), 
C(u 5 ), S(u 2 ), S(u 5 ) terms in the first factor of Eq. (11) are non- 
zero, with u 2 =3a/2 and w 5 =-3a/2. In this case, Eq. (11) reduces 
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to Eq (6) for a single aperture with an aperture width of 3a, as 
it must. 

From Eq.(lO) or Eq.(ll), it is clear that the calculation of 
electric field or intensity at a point P requires, in general, the 
evaluation of 16 Fresnel integrals. 12 of these involve the 
argument in u (or y) only, and the remaining 4 in v (or z) only. 
These equations are the basis of calculation of the complete 
intensity distribution of the Fresnel diffraction pattern from a 
triple aperture, as will be seen in the next section. Either of 
these two equations can be used, but in this work Eq. (10) was 
chosen to calculate firstly the electric field, and then the 
intensity was calculated by taking the square of it. These 
equations have two factors, the first factor involves only the u 
(or y) coordinate and expresses the dependence of electric 
field or intensity in the y direction. The second factor involves 
only the v (or z) coordinate and expresses the dependence of 
electric field or intensity in the z direction. 




Fig. 2 Configuration of the triple -aperture problem for Fresnel diffraction, a 
is the individual aperture width and (a+b) is the aperture separation. 

ni. Simulation technique 

Eq. (10) and Eq. (11) describe the electric field and the 
intensity at P, respectively. In order to describe these for an 
off axis point P', a method outlined in [4] and [5] was used. 
The observation screen and SOP line were fixed. Then instead 
of moving P, the aperture in yz plane was moved in opposite 
direction, so that the relative position of the aperture is this 
new position and point P remains unchanged. We then 
calculate the electric field at P instead of at P'. Point P' will 
see a new set of values for y's and z's (and therefore for u's 
and v's). For example, to find the intensity at point P' 1mm to 
the left of P, the screen was kept undisturbed and the aperture 
was moved 1mm rightwards, and the intensity at P in this 
configuration is calculated instead of at P'. Consequently, the 
point P will now see a new set of values for y h y 2 , y3, y^ ys 
and y 6 and therefore, for Uj, Il2, U3, U4, U5 and U6 in Eq. (10) 
and Eq. (1 1). In principle, the electric field and the intensity at 
any point P' on the image plane can be found in this way by 
making appropriate translations of the aperture in the y and z 
directions, and in Eq. (10) or Eq. (12), substituting 
correspondingly a new set of values of u and v limits. 

The flow chart of the algorithm is shown in Fig. 3. For 
calculation of the electric field distribution using Eq. (10) at 
all the points (pixels) on the image plane, a large number of 
the Fresnel cosine and sine integrals will need to be evaluated 



efficiently. This is accomplished in the MATLAB program 
using the special functions mfun ('FresnelC, i:s:f ) and 
mfun ('FresnelS', i:s:f ). For example, the MATLAB 
statement R=mfun ('FresnelC', i:s:f ) generates an array R of 



Anpi 



Input parameters amm, bmm, cram 
Witim, smm, qOmm, Inm, 1 in a GUI 



Convert dimensional variables amm, 
bmm, cmm, Wmm, smm into 
dirnenionless variables, a,b,c,W,s 



X 



For the first quadrant, define limits of 
Fresnel integrals 3a/2+b to W+3a/2+b 
for u 2j a/2+b to W+a/2+b for u n , a/2 to 
W+a/2 for u 4 , -a/2 to W-a/2 for u 3 , -a/2-b 
to W-a/2 -b for u 6 and -3a/2-b to 
W-3a/2-b for u 5s and calculate arrays 
Fresnel integrals C(u 2 ), C^), C(u 4 ) 
C(u 3 ), C(u fi ), C(u 5 ), S(u 2 ), S( Ul ), 
S(u 4 ), S(u,), S(u 6 ), S(u 5 ), S(u 2 ), S( Ul ) 
for y (or u) axis 



Calculate 1 -D array A indicating 
variation of complex E-field in y(or u) 
direction 



For the first quadrant, define limits of 
Fresnel integrals b to W+b for v a and-b 
to W-b for V, and calculate arrays of 
Fresnel integrals Qv^Qv^SfV;, ),S(v,) 
for z (or v)-axis 



Calculate l-D array B indicating 
variation of E-field in z (or v)-axis 



Construct square matrix B whose 
element values vary in row direction 



1 



Construct square matrix A whose 
element values vary in the column 
direction 



Multiply B & A to construct 2-D 
matrix C to indicate E-field 
changes in u or v axes for first quadrant 



Square matrix C to calculate the intensity 
variation in both axes for first quadrant 



Normalize all elements value of matrix 
D and multiply by intensity factor t 



Construct 4 -quadrant zero matrix E 



Place elements of D into first 
quadrant of E 



Invert E around z-axis 



Invert E around y-axis to construct 
a 4-quadrant image matrix 



Construct l-D array y indicating 
real position in mm 



Show 4-quadrant intensity 
as an image in real dimensions 
hnagesc command 



matrix E\ 
ions by J 



Fig. 3 Flow chart of the algorithm for the triple aperture problem 

the Fresnel cosine integrals with arguments staring from i and 
ending in /, with a step size s. Since the triple aperture is 
symmetrical about the y and z axes, the Fresnel diffraction 
pattern generated by it will obviously be symmetrical about 
the Y and Z axes in the image plane. Therefore, by using this 
inversion symmetry around Y and Z axes, the amount of 
calculation can be reduced to one-fourth. Therefore, the 
diffraction image in only one quadrant (e.g. the first) in the 
YZ plane needs to be generated, and then the image in the 
second quadrant can generated by simply inverting it around 
the Z-axis. Another inversion around the Y-axis of these two 
image quadrants will produce the images in the third and 
fourth quadrants. 

As shown in Fig. 4, the aperture was displaced (virtually) 
by an amount W, and the extreme limits of the edges of the 
displaced aperture were determined and compared to the 
corresponding limits for the undisplaced aperture to find the 
range of u and v values. When the aperture was moved by W 
to the left, the ranges for the six edges of the triple aperture 
will be established as: (3a/2+b) to (W+3a/2+b) for y 2 , 
(a/2+b) to ( W+a/2+b) for y h (a/2) to (W+a/2) for y 4 , (-a/2) to 
(W-a/2) for y 3 , (-a/2-b) to (W-a/2-b) for y 6 and (-3a/2-b) to 
(W-3a/2-b) for y 5 . Arrays of both the Fresnel cosine and sine 
integrals need to be calculated corresponding to these six 
ranges by the mfun statements. By using these 12 arrays, 
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the electric field or intensity dependence in the y (u) 
direction can be calculated (see Eq. lOorEq. 11). Following 



z 




Y 



Fig. 4 Apparent limits of the displaced triple aperture seen from the 
observation area. 

a similar procedure for the entire aperture, i.e. by moving the 
aperture by W downwards, the zi and z 2 ranges are determined: 
c to (W +c) for Z2 and -c to (W-c) for zi. Four arrays of 
Fresnel integrals are to be evaluated for these two ranges, 
providing numerical values for the calculation of the second 
factor in Eq. (10). 

After the initial parameters to the program are input 
through a GUI (Graphical User Interface) generated in lines 1- 
2 of the program, the dimensionless quantities corresponding 
to aperture dimensions a, b and c, step size s and image area 
size W are calculated in lines (4-5) of the program. The step 
size s determines the resolution of the simulated images, since 
the number of pixels in the image is given by 2W/sX2W/s. 
The 12 Fresnel integrals arrays for the u (or y) dimension are 
evaluated in lines (6-17). In line 18 A the electric field in y (u) 
is calculated in complex form, corresponding to the first 
complex factor in Eq. (10). In the next four lines, the Fresnel 
cosine and sine integrals are evaluated for the z (or v) 
dimension. Next, in line 23, the second complex electric field 
in v is calculated, corresponding to the second factor in Eq. 
(10). In the next lines, the matrix C, which contains the 
complex electric field variation in both u and v directions, is 
constructed. Finally, by squaring it, the matrix D, which 
contains the intensity variations for the first quadrant of the 
image plane, in both y and z, is calculated (line 26). The 
elements of this matrix then normalized (line 27), and the 
complete E matrix for the full image plane intensity 
distribution is constructed by inverting D twice (lines 30-38). 
Finally, the image is displayed as a greyscale image by the 
images c command (line 40). The complete MATLAB 
program, which we call tpslit, is given in the Appendix. 

IV. Simulation results 

Using the program tpslit, the values of the inputs, i.e. 
aperture parameters a and b in mm (see Fig. 2), height c in mm, 
image area W in mm, step size s in mm, the aperture image- 



image plane distance q in mm, the illuminating wavelength / 
in nm and the intensity factor t are input by a GUI to the 
program (Fig. 5a). The intensity factor t is used to control the 
apparent visual intensity in the generated image. The 
computer generated Fresnel image due to a triple aperture for 
wavelength X=500 nm, q =400 mm, image area W=4 mm, 
aperture width <2=0.5 mm, b=2 mm (with aperture separation 
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(a) 




-4-3-2-101234 



(b) 

Fig. 5(a) Screenshot of the MATLAB GUI, showing the 8 input parameters 
(b) The computer simulated Fresnel image from triple aperture. 

2.5mm), height c=3mm and step size s=0.0 1mm is shown in 
Fig. 5b. Three separate and distinct Fresnel diffraction 
patterns can be observed from each aperture, separated by a 
distance, with little mutual interference of the diffracted light 
between the three individual apertures. 
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(a) (a+b)=2.5mm 



(b) (a+b)= 1.5mm 





(c) (a+b)= 1.0mm 



(d) (a+b)=0.5mm 



Fig. 6 Computer simulated Fresnel diffraction images for decreasing (center-to-center) aperture separation (a+b), while keeping the aperture width (a=0.5 mm) 
constant. The wavelength is 500 nm, the aperture-screen distance q is 400 mm, step size s=0.01 mm and c=3 mm 





(a) a=0.3mm 



(b) a=0.2mm 





(c) a=0.1mm 



(d) a=0.05mm 



Fig. 7 Computer simulated Fresnel diffraction images for decreasing aperture width a, while keeping the aperture separation (a+b) constant at 1.5 mm. The 
wavelength is 500 nm, the aperture -screen distance is 400 mm, step size s=0.01 mm and c=3 mm. 
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Next, in order to examine the effect of a change of 
simulation parameters on the diffracted images, we have 
generated a series of images. In the first set of simulations, the 
aperture separation (a+b) was made successively smaller, 
while keeping the aperture width (a=0.5mm) constant. A 
series of diffraction images were simulated as shown in Fig. 
6(a)-(d). 

In Fig. 6(a), in the region between the apertures, there was 
not much interference between the light from each of the 
aperture diffraction patterns. But, in Fig. 6(b) and (c), some 
distinct interference between diffracted light from the three 
individual apertures has been observed as the separation of the 
apertures is reduced. As the separation was decreased to zero 
[fig. 6(d) for b=0], a typical single aperture Fresnel diffraction 
pattern from an aperture of dimensions 1.5mmX3mm has 
been observed. 

In the next series of simulations, the width a of the 
apertures was made gradually narrower while their separation 
(a+b) was kept constant at 1.5 mm. This is shown in Fig. 
7(a)-(d). From Fig. 7(a) and (b), as the apertures become 
narrower, the light is diffracted over a wider region and some 
interferences between diffracted light occur. In Fig. 7(c) and 
(d), at smaller aperture widths, mutual interference between 
diffracted light produced fringes which look like typical 




Young's fringes. Light from each aperture is diffracted over a 
wide portion of the image area and interfere with each other. 
In Fig. 7(c) and (d), the density of the 'fringes' is the same, 
since this depends only on aperture separation. 

V. COMPARISON WITH FRAUNHOFER DIFFRACTION 

If the aperture-screen distance q is increased, a gradual 
transition from the Fresnel regime to Fraunhofer regime 
should be expected. As a rule of thumb, if D is the dimension 
of the square aperture, then Fresnel diffraction occurs if D 
satisfies the following relation [2,10] , 

q < (D 2 /A). (12) 

Otherwise, Fraunhofer diffraction should occur. In the 
next series of simulations, we increased aperture-screen 
distance q , while keeping all other parameters constant [Fig. 
8]. A transition from the Fresnel regime to Fraunhofer regime 
is observed, being consistent with Eq.(12). For example, in 
Fig. 8(a), if we take the lateral dimension to be D=lmm, then 
D /X is about 2000 mm, and we can expect the diffraction 
to be Fresnel-like in the ^-direction. On the other hand, in Fig. 
8(d), (D 2 /X) >q , and we can expect the diffraction pattern to 
be Fraunhofer-like. In between these extremes, in Fig. 8(b) 
and 8(c), a transition from Fresnel to Fraunhofer regime is 
seen to occur. 




-5-4-3-2-1 1 2 3 4 5 



(a) q =400mm 



(b) q =800mm 




(c) q = 1500mm (d) q =8000mm 

Fig. 8 Computer simulated Fresnel diffraction images for increasing aperture-screen distance q , while keeping all other parameters constant, with a=0.5 mm, 
b=0.5 mm, c=3 mm. The wavelength is 500 nm and step size s=0.01 mm. (a) is in the Fresnel regime while (d) is in the Fraunhofer regime, according to Eq.(12) 
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(a) )i=500nm 



in 



(b) X=1000nm 



imi 



(c) ^=2000nm 



(d) ^=10000nm 



Fig. 9 Computer simulated Fresnel diffraction images for increasing wavelength X, while keeping all other parameters constant, with a=0.5 mm, b=0.5 mm, 
c=3 mm and q =400 mm. (a) is in the Fresnel regime while (d) is in the Fraunhofer regime, according to Eq.(12) 



By using the program tpsllt, it is also possible to 
observe the effect of change of wavelength on the diffraction 
patterns. For example, when the wavelength was increased 
from 500nm (green in the visible region) to lOOOOnm, while 
keeping everything else constant in the experimental 
configurations (including q ), a transition from the Fresnel 
regime to Fraunhofer regime should occur, according to Eq. 
(12). This is shown in Fig. 9. In Fig. 9(a), the diffraction 
pattern is definitely Fresnel-like, while that in Fig. 9(d) is 
Fraunhofer-like. At intermediate wavelengths (e.g. at 
A,=1000nm and A,=2000nm), a transition from Fresnel to 
Fraunhofer regime is actually observed. 

The simulated Fraunhofer diffraction image for 
g =8000mm, and an aperture dimensions a=0.5mm, b=0.5 
mm and c=lmm, and /l=500nm is shown in Fig. 10(a). For 
quantitative comparisons, the three dimensional mesh plot of 
this image is shown in Fig. 10(b), which is generated by the 
MATLAB mesh command. In this situation, the value of D 2 /X 
is 2000mm, and we expect a transition to Fraunhofer regime 
to occur since q is more than this value. 

The Fraunhofer diffraction pattern of a multiple aperture 
in the image plane (Y, Z) can be calculated by an analytical 
formula [1,2] 



/ ( r,Z) = 7 (^)(^)(^). (13) 



sin y 



a 



where N is the number of apertures (Af=3 in our case) and ft, y 
and a are defined as, 



R n Y 
B = —a—. 
X R 



Y = — (a + b) — . 
' X R 



n Z 
a = —c — . 
X R 



(14) 



In the above equations, a is aperture width, (a+b) is the 
aperture separation, c is aperture height, X is the wavelength 
and R is the distance between the aperture and the screen, 
assumed to be very large. Using Eq. (13), a MATLAB 
program was written to compute the intensity distribution in 
the image plane, which accepts the values of a, b, c, distance 
R, wavelength X and the desired image area 2WX2W. The 
intensity can be shown, for quantitative comparisons, as a 
mesh plot by MATLAB. The mesh plot is shown in Fig. 10(c) 
considering the above parameters. 

Comparing Fig. 10(b) and 10(c), a fair agreement can be 
observed, including the minor peaks. This means that the 
simulation method of Fresnel diffraction for triple aperture is 
correct, since it has a good quantitative match with a purely 
Fraunhofer calculation under similar conditions, i.e. at large 
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aperture- screen distance, q . The slight difference between Fig. 
10(b) and 10(c) may be attributed to the fact that Eq. (13) is 
valid for extremely large distances, but in reality the distance 
is 8000mm only. 

To test this assumption, we have generated another 
computer simulation by our program for the above situation, 
but with a larger value of q (20,000mm). The intensity 
distribution and mesh plot are shown in Fig. 11(a) and 11(b). 



Fig. 11(c) shows the mesh plot of the intensity distribution 
calculated by the Fraunhofer formula [Eq. 13] for g =20,000 
mm. In this case, a better agreement between Fig. 11(b) and 
Fig. 11(c) can be clearly seen. This means that as q becomes 
larger and larger, the agreement between the simulation 
results by IFIM method and the Fraunhofer formula becomes 
better and better. 





(c) 

Fig. 10 (a) Computer simulated far-field diffraction images for the triple aperture, with a=0.5 mm, b=0.5 mm, c=l mm, X=500 nm and q = 8000 mm (Fraunhofer 
limit), (b) Mesh plot of the far-field diffraction pattern for this aperture, (c) Mesh plot of the far -field diffraction pattern for this aperture generated by the 

Fraunhofer formula. 
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(c) 

Fig. 11 (a) Computer simulated far-field diffraction images for the triple aperture, with a=0.5 mm, b=0.5 mm and c=l mm, X=500 nm and q = 20,000 mm 
(Fraunhofer limit), (b) Mesh plot of the far-field diffraction pattern for this aperture. (c) Mesh plot of the far-field diffraction pattern generated by the 

Fraunhofer formula. 



VI. DISCUSSION 

In this paper, we have described a method to calculate and 
simulate the complete Fresnel diffraction image of the triple 
aperture using the Iterative Fresnel Integrals Method. Though 
not as popular as the double aperture or slit, the triple slit has 
been used as the source of interference in a triple-slit 
interferometer [11]. Each of our simulations is a virtual 
experiment, where the results of performing a real diffraction 
experiment can be obtained purely by simulation, without 
actually performing it. 

The theoretical background as well as the complete 
implementation of the algorithm in MATLAB codes is given 
in this paper. We have computed first the complex electric 
field distribution in the image plane, and then calculated the 
intensity distribution by squaring it. Here we have used the 
complex number capability of MATLAB to full advantage. 
This affords some simplification in the computation process. 

In Eq. (2) or Eq. (3), it was tacitly assumed that the 
obliquity factor K(<f>) is unity for all waves emitted from all 
the surface elements dS located inside the triple aperture. But 
actually, in our simulations, in order to calculate the electric 
fields at an arbitrary point P' in the image plane, we virtually 
translated the triple aperture in the y and z directions and 
calculated the electric field at the origin P of the image plane 
(see, for example, Figs. 2 and 4). These translations of the 
aperture effectively make the light emitted from the aperture 
slightly non-paraxial (off-axis). As an example, [in the 
simulation of Fig. 5], for an image area W=4 mm, aperture- 
screen distance q =400mm, and with aperture dimensions 
£/=(). 5 in in, Z?=2mm, c=3mm, the extreme off-axis waves will 
be emitted from surface elements dS located at y=W+3a/2+b 
and at z=W+c, for virtual translations Wof the aperture in the 
y and z-directions, respectively (see, for example, Fig. 4). For 
these waves, the maximum values of (/> are (W+3a/2+b)/q 
and ( W+c)/q , respectively, and work out to be about 1 degree 
in both cases. The corresponding value of the obliquity factor 
K((f>) [=V6(l+cos (|>)] is about 0.99992, and can be taken equal 
to unity, as was done in Eq. (2) and Eq. (3). This 
approximation produces a maximum error of less than 0.01% 
in the calculation of electric fields. For larger image areas, 
larger apertures, or smaller aperture-screen distances, (/> would 



be correspondingly larger. To limit the error in K(<pJ (and 
hence in the electric field) to less than 1%, the maximum 
allowed value of ^should be about 10 degrees. This probably 
covers most cases of practical interest. 

For the simulations, we used an office computer using an 
Intel core-i7 CPU having 4GB of RAM running at a clock 
speed of 3.4GHz under Windows 7. The simulation time for a 
typical 1000X1000 pixel image is about 10 seconds. For 
slower computers or for larger images, the simulation time 
may be significantly greater. 

Using the program, we simulated the Fresnel diffraction in 
a number of situations. We have found that as the individual 
apertures of the triple aperture are brought closer, significant 
interference between diffracted light can be seen in the region 
between the apertures. On the other hand, if the aperture 
widths of the individual aperture are made narrower while 
keeping their separation constant, then the light is diffracted 
over a wider region. Interference effects with the light 
diffracted from the apertures can be seen, and Young-like 
fringes appear on the image plane. 

The expected transition from Fresnel to Fraunhofer 
regions can also be observed under appropriate conditions. 
For example, when the aperture-screen distance q is 
increased, or when the wavelength X is increased, a transition 
from Fresnel regime to Fraunhofer regime is actually 
observed. In addition, good quantitative agreement was 
observed between the simulated Fresnel diffraction pattern at 
large aperture-screen distances (in the Fraunhofer regime) and 
the intensity distribution generated by an analytical 
Fraunhofer formulation. This agreement becomes better as the 
aperture-screen distance is increased more, and provides a 
proof that our simulation method is correct and appropriate. 

VII. CONCLUSION 

In conclusion, we have described a method to compute 
and generate, in virtually any PC, the complete near-field 
Fresnel diffraction patter from a triple aperture in any 
arbitrary experimental configuration by the Iterative Fresnel 
Integrals Method. This simulation and virtual experiment 
provide a simple and easy-to-use method to study the complex 
phenomenon of diffraction from a triple aperture system. 
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APPENDIX 

tpslit: Complete MATL AB program for the triple 
aperture 

1. u=inputdlg({'a mmVb mmVc mmVW mmVs mm',' 

qO irimyi nm',... 

2. 'exposure' },'Fresnel Diffraction from triple aperture', 

[1,1,1,1,1,1,1,1]); 

3. for i=l:8; v(i)=str2num(u{i}); end 

4. q0=v(6); t=v(8); l=v(7)*le-6; f=sqrt(2/(l*q0)); 

5. a=v(l)*f; b=v(2)*f; c=v(3)*f; W=v(4)*f; s=v(5)*f; 

6. Cu2=mfun('FresnelC',3/2*a+b:s:W+3/2*a+b); 

7. Cul=mfun('FresnelC',a/2+b:s:W+a/2+b); 

8. Cu4=mfun('FresnelC',a/2:s:W+a/2); 

9. Cu3=mfun('FresnelC',-a/2:s:W-a/2); 

10. Cu6=mfun('FresnelC',-a/2-b:s:W-a/2-b); 

11. Cu5=mfun('FresnelC',-3/2*a-b:s:W-3/2*a-b); 

12. Su2=mfun('FresnelS',3/2*a+b:s:W+3/2*a+b); 

13. Sul=mfun('FresnelS',a/2+b:s:W+a/2+b); 

14. Su4=mfun('FresnelS',a/2:s:W+a/2); 

15. Su3=mfun('FresnelS',-a/2:s:W-a/2); 

16. Su6=mfun('FresnelS',-a/2-b:s:W-a/2-b); 

17. Su5=mfun('FresnelS',-3/2*a-b:s:W-3/2*a-b); 

1 8 . A=complex(Cu2+Cu4+Cu6-Cu 1 -Cu3 - 

Cu5,Su2+Su4+Su6-Sul-Su3-Su5); 

19. Cv2=mfun('FresnelC',c:s:W+c); 

20. Cvl=mfun('FresnelC',-c:s:W-c); 

21. Sv2=mfun('FresnelS',c:s:W+c); 

22. Svl=mfun('FresnelS',-c:s:W-c); 

23. B=complex(Cvl-Cv2,Svl-Sv2); 

24. B=B';n=size(B);B=repmat(B(:,l),l,n); 

25. A=A';A=repmat(A(:,l),l,n);A=(A)'; 

26. C=B*A;D=C.*conj(C); 



27. D=t*D/max(max(D)); 

28. m=2*fix(((2*W)/s)/2); 

29. E=zeros(m+l,m+l); 

30. for q=l:l:m/2+l; 

31. for p=l:l:m/2+l; E(m/2+2-p,q+m/2)=D(p,q); end 

32. end 

33. for q=m+l:-l:m/2+2; 

34. for p=l : 1 :m/2+l ;E(p,-q+2+m)=E(p,q);end 

35. end 

36. for q=l:l:m+l; 

37. for p=l:l:m/2;E(-p+2+m,q)=E(p,q);end 

38. end 

49. y=-W:s:W;ymm=y/f; 

40. imagesc(ymm,ymm,E,[0,l]);colormap(gray); 
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